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Thermal convection in nanofluids is investigated by means of a continuum model for binary-fluid 
mixtures, with a thermal conductivity depending on the local concentration of colloidal particles. 
The applied temperature difference between the upper and the lower boundary leads via the Soret 
effect to a variation of the colloid concentration and therefore to a spatially varying heat conductivity. 
An increasing difference between the heat conductivity of the mixture near the colder and the warmer 
boundary results in a shift of the onset of convection to higher values of the Rayleigh number for 
positive values of the separation ratio i[> > and to smaller values in the range ip < 0. Beyond some 
critical difference of the thermal conductivity between the two boundaries, we find an oscillatory 
onset of convection not only for tp < 0, but also within a finite range of ij) > 0. This range can 
be extended by increasing the difference in the thermal conductivity and it is bounded by two 
codimension-2 bifurcations. 

PACS numbers: 47.55. P-, 47.57.E-, 44.10. +i 
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I. INTRODUCTION 

Thermal convection plays a central role in geophysics 
[H> > atmospheric dynamics Q and various technical ap- 
plications, whereof nanofluids were recently identified as 
efficient heat-transfer substances 0-0|- Rayleigh-Benard 
convection with its numerous variants is also a classi- 
cal lab experiment for studying generic phenomena of 
nonlinear dynamics and pattern formation, which occur 
in this system as in several other fields of natural sci- 
ence @, 18-flfl] . In most of the related theoretical studies 
the Oberbeck-Boussinesq (OB) approximation for ther- 
mal convection is used, where constant material parame- 
ters independent of the thermodynamic variables are as- 
sumed, except the temperature-dependent density within 
the buoyancy term, which is the essential driving force of 
convection. 

Non-Boussincsq contributions to the governing equa- 
tions for convective systems are often required to model 
various phenomena in an appropriate way. Around 4°C, 
for instance, the linear term in the thermal expansion 
of water vanishes and the quadratic contribution has to 
be taken into account. It is well known that this mod- 
ification changes the symmetry condition in a thin fluid 
layer heated from below, leading to hexagonal convection 
patterns instead of stripe patterns close to the onset of 
convection [H, [ll[ . Strongly varying material properties 
in the Earth's mantle arc a major motivation for using 
a temperature-dependent viscosity in models of thermal 
convection in single component fluids, which is another 
non-Boussinesq contribution [12h17| . A further example 
is a temperature dependent thermal conductivity, being 
considered to be important, for instance, for explaining 
a delayed cooling of the Earth's mantle There are 
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also recent studies about non-Boussinesq contributions to 
convection in binary-fluid mixtures where either a tem- 
perature dependent thermo-diffusion coefficient or a de- 
pendence of the viscosity on the local composition was 
taken into account [19l - [2ll |. 

Convection in binary-fluid mixtures is another classi- 
cal, driven pattern- forming system [22|, [23[, which has 
attracted wide attention during the last decades (22l - [30| . 
In binary-fluid mixtures the concentration field of one of 
the two constituents enters the basic equations as an ad- 
ditional dynamic quantity [22T - I241 ] . Via the Soret effect 
(thermophoresis) a temperature gradient, that is applied 
vertically across a convection cell, may cause variations 
of the concentration field which couples into the Navier- 
Stokes equations for the velocity field via the buoyancy 
term. Depending on the sign of the Soret effect, the 
heavier constituent is either driven to the colder upper 
boundary or to the warmer lower boundary. In the for- 
mer case one obtains stationary convection patterns near 
the onset of convection and in the latter case oscilla- 
tory patterns. Experimentally, the onset of convection 
is well investigated in mixtures of alcohol and water as 
well as for 3 He/ 4 He mixtures Q- The possibility of hav- 
ing both, a stationary as well as an oscillatory onset of 
convection, including a so called codimension-2 bifurca- 
tion at the transition point, made it to a very attractive 
model system for generic bifurcation phenomena [g, |31| . 

Colloidal suspensions, also known as nanofluids, may 
be considered as a further example of a binary mix- 
ture, with the suspended particles being the second con- 
stituent. Recently, convection in colloidal suspensions 
[32T - |36| has been investigated experimentally with a spe- 
cial focus on Soret driven convection [32h34| | , on bistable 
heat transfer which is caused by sedimentation effects 
i35] , or on the effects of thermosensitive particles (3(| . 

Additionally, in several nanofluids with particle sizes in 
the range of 1 — 100 nm a strong dependence of the ther- 
mal conductivity on the concentration of nanoparticles 
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was reported, see, e.g., Refs. d @, 0, UMl- This de- 
pendence presents another non-Boussinesq effect. In this 
work we investigate its impact on the onset of convec- 
tion. Taking the particle-concentration-dcpcndcnt ther- 
mal conductivity into account, we go beyond other recent 
studies with the focus on convection in colloidal suspen- 
sions [40l l4l| , where the Boussinesq approximation has 
been used. 

The enhancement of the thermal conductivity in fluids 
by increasing the concentration of suspended nanopar- 
ticles was confirmed by a recent benchmark study with 
contributions of more than thirty researchers [6| . In spite 
of the promising applications of nanofluids for improving 
heat transfer in cooling systems, a considerable number 
of open questions is left. Among several possible trans- 
port phenomena, discussed to explain the heat transfer 
enhancement in nanofluids, Brownian diffusion and ther- 
mophorcsis were identified as the two most important 
ones [5j. Both mechanisms build the crucial extensions 
from models for a single component Newtonian fluid to 
models for convection in binary fluid mixtures [HI, [24| . 

The paper is organized as follows: In Sec.|H]we briefly 
present the underlying equations of motion and introduce 
a linear relation between the heat conductivity and the 
particle concentration, which leads to a nonlinear spa- 
tial dependence of the heat conductive state. The es- 
sential methods of the linear stability analysis for the 
determination of the onset of convection are presented 
in Sec. Mil Our numerical results for the onset of con- 
vection include a prediction of an oscillatory onset even 
in the range of positive values of the separation ratio for 
both cases, no-slip boundary conditions in Sec. MI Al and 
free-slip boundary conditions in Sec. lIIIBl In Sec. II VI we 
summarize and discuss our results. 



II. BASIC EQUATIONS AND HEAT 
CONDUCTIVE STATE 

To describe convection in a horizontal layer of a col- 
loidal suspension the common mean field approach for 
binary-fluid mixtures is used [H, [H H, [H, [Hi El] . In 
addition, we take into account a linear dependence of 
the thermal conductivity k on the mass fraction of the 
colloidal particles N(r,t). The mass density of the col- 
loidal particles p c is assumed to be similar to the mass 
density of the solvent p s , i. e. e = p c /p s — 1 • Fur- 
thermore, we assume small colloidal particles undergoing 
strong Brownian motion, so that sedimentation effects 
become negligible. Deviations of N(r, t) from the mean 
mass fraction No may lead to a spatial dependence of the 
thermal conductivity via the linear relation 



k (1 + 7 (N-N Q )) , 



(1) 



where kq describes the mean thermal conductivity of the 
suspension and 7 = (1/kq) dn/dN is a measure for the 
dependence of the thermal conductivity on the concentra- 
tion of the nanoparticles. Heat conduction experiments 



with small volume fractions of nanoparticles match with 
the assumption 7 ~ 2.5 d [H E^. 

The common set of basic transport equations for in- 
compressible binary- fluid mixtures, cf. Refs. [H, HE HH, 
HslHilj involves the temperature field T(r, t), the mass 
fraction of the particles N(r, t), the fluid velocity v(r, t), 
the density of the mixture p(r, t), and the pressure field 
p(r,t): 



V • v = 0, 
{d t + v • V) T = V • ( X VT) , 

(d t + v • V) JV = D V • ( VN + VTJ , 



{dt + v • V) v 
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(2a) 
(2b) 

(2c) 
(2d) 



Eq. ([2a)) describes the incompressibility of the fluid. \ = 
k/ po in the heat equation (|2bl) denotes the thermal diffu- 
sivity of the mixture, which is in our model a function of 
the concentration of the suspended colloids. D in Eq. (|2c|) 
is the diffusion constant, which takes due to the size of 
the colloidal particles much smaller values than in molec- 
ular binary-fluid mixtures. The dimensionless thermal- 
diffusion ratio fc<r representing the cross coupling between 
the temperature gradient and the particle flux is related 
to the Soret coefficient S T via k T /T = N(l - N)S T , 
which can be either positive or negative. Throughout 
this work fcp/T ~ Nq(1 — Nq)St is regarded as constant. 
v in the Navier-Stokes-equations (|2dl) is the kinematic 
viscosity. The gravity field g = — ge z is chosen parallel 
to the z-direction. 

For the local density p of the suspension we use a lin- 
earized equation of state [13, SH , 



p = po [1 - a(T - To) + I3(N - N )} , 



(3) 



with the thermal expansion coefficient a = 
-{l/pa)dp/dT and j3 = (l/p )dp/dN reflecting the 
density contrast between the solvent and the suspended 
particles. According to the Boussinesq approximation 
this dependence of the density is taken into account 
only within the buoyancy term. The sign of /3 indicates 
whether the colloidal particles have a higher or a lower 
mass density compared to the solvent. Here, we assume 
P > corresponding to e > 1 . 

Boundary conditions. The fluid layer is confined be- 
tween two impermeable, parallel plates at a distance d, 
and extends infinitely in the (a;-y)-plane. The lower plate 
at z — —d/2 is kept at a higher temperature To + 6T/2 
than the upper plate at z = +d/2 with the lower temper- 
ature To — 5T/2. Together with a vanishing mass current 
at the boundaries and realistic no-slip conditions for the 
flow field we have the following set of boundary condi- 
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tions at z = ±d/2: 



T = T oT 



ST 
2 ' 



= d z N+-^d z T, 
Jo 

= v x = v y = v z = d z v z 



(4a) 

(4b) 
(4c) 



For geophysical applications free-slip boundary condi- 
tions 



= d z v x = d z v y = v z = d z v z , 



(5) 



for the flow field are often considered to be more realis- 
tic [HI. In the case of a constant thermal conductivity 
and free-slip, permeable boundary conditions an analyt- 
ical determination of the onset of convection is possible 
25] . By introducing a concentration dependent thermal 
conductivity this advantage is lost and one has to rely on 
numerical methods. 



A. Heat conductive state 

Together with the concentration dependent heat con- 
ductivity given by Eq. ([1]) the constant vertical heat 
current j z = — x<9 z T con( j leads to a nonlinear z- 
dependence of the time-independent temperature distri- 
bution T con( j(z) and the corresponding concentration dis- 
tribution N ccm d(z) of the heat conductive state. Both are 
derived in the Appendix „ and are of the form 

T cond = T + ^ Q (1 + Y(0) - W(z, (6a) 
= T -^-4(i-JW), (6b) 



N cond = N + ^-(l-W(z,0 



(6c) 



with the abbreviations 



w(z, o = M (i + Y(0) 2 + ~ (i + HO + ^ 2 ) • 



and 



SN = -^ST. 
For finite values of 7 the vertical heat current, 
XoST 



3z 



2d 



1 



(7a) 
(7b) 

(7c) 
(8) 



is reduced compared to the situation of a constant heat 
conductivity being independent of the particle concen- 
tration. 

We would like to stress that Eqs. © are derived under 
the assumption of small values of SN, but these formulas 
are still reasonable for values up to SN < No/3. The re- 
striction |£| < v3 according to Eq. (|7bl) has the same ori- 
gin and is fulfilled by all parameters chosen in this work. 
For stronger variations of SN the assumption of a con- 
stant thermodiffusion coefficient kr/T ~ Nq(1 — No) St 
is not justified anymore and a generalized approach has 
to be chosen. 



B. Dimensionless equations of convective fluid 
motion 

For the further analysis it is convenient to separate 
the basic heat conductive state in Eq. ([5]) from convec- 
tive contributions as follows: T(r, t) — T cond (z) +T%(r, t) 
and N(r,t) = N cond (z) + Ni(r,t). Using the rotational 
symmetry in the fluid layer we can restrict our analysis to 
two spatial dimensions, namely to the (x-z) -plane. With 
this simplification the fluid velocity v = (v x ,0,v z ) can 
be expressed by a stream function <f>(x, z, t): 



9.r 



-d z 



(9) 



Subsequently all lengths are scaled by the vertical dis- 
tance d and times by the vertical thermal diffusion time 
d 2 /xo- Scaling the temperature field T by (xo^o) /(&gd 3 ), 
the concentration field N by — (hrXoVo) / (Toagd 3 ) and 
the stream function </> by Xod we are left with five dimen- 
sionless parameters: The Rayleigh number R, the Prandtl 
number P, the Lewis number L, and the separation ratio 

4>, 



P= V JL. L = R. R 



Xo 



Xo 



a ^. ST , tf=^, (10) 



Xo"o 



aT 



are well known from molecular binary-fluid mixtures \2t 
[29|]. The fifth dimensionless quantity, 



C = 7 



^0X0 



(11) 



is introduced to characterize the spatially varying con- 
tribution to the thermal diffusivity caused by the 
concentration-dependence of the thermal conductivity of 
the suspension. An illustration of its physical meaning 
is obtained by considering the thermal conductivity con- 
trast between the upper and the lower boundary 



l + ^ 7 fc T ST /Tp _ l + \R^C, 
1- ±jk T ST/T ~ 1 -hRil>C 



(12) 



which is essentially a function of the product of the three 
dimensionless control parameters Rip( = £. 

In the following we will discuss our results essentially in 
dependence on tp and £, whereas P and L are considered 
as constants. 
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Introducing a rescaled temperature deviation 9 — 
(R/ST)Ti, a rescaled concentration deviation Ni — 
— (T Q R/kTST)Ni, and a rescaled stream function $ — 
l/(Xod) 4> m terms of these dimensionless quantities and 
using the combined function c = N\ — 9 instead of N\ we 
obtain: 

d t -W A9-Rip (Wi(d z c+2d z 9) (13a) 

- W 2 (c + 9) - RWid x § = 

- i>( {(d x 8)d x {6 + c) + (d z 9)d z (9 + c) + A9{9 + c)] 
+ {d z $d x -d x W z )8, 

d t c + W A9 + R%PCWi{d z c + 2d z 9) (13b) 
+ W 2 (c + 9) - LAc = 

+ tP( \{d x 9)d x {9 + c) + {d z 9)d z {9 + c) + A9{9 + c)] 
+ {d z $d x -d x W z )c, 

d t A<S> - PA 2 $ - P(l + ^)d x 9 - Ptpd x c= (13c) 
P(d z $d x - d x $d z ) A$ , 

with the abbreviations 

W = W(z,0, (14a) 

^1=1^,0, (14b) 

W 2 = dlW{z,t). (14c) 

For reasons of simplicity all the tildes have been sup- 
pressed. No-slip, impermeable boundary conditions for 
the fields 9, c, and 5> demand 

9 = d z c = $ = <9 2 $ = at z = ±~ (15) 

while free-slip, permeable boundary conditions [III |42| 

= c = $ = 9f $ = at z = ±- . (16) 

In this work we present both, results based upon no- 
slip, impermeable boundary conditions and results based 
upon free-slip, permeable boundary conditions. 

III. LINEAR STABILITY OF THE HEAT 
CONDUCTIVE STATE AND ONSET OF 
CONVECTION 

In order to determine the parameters at the onset of 
convection we investigate the linear stability of the heat 
conductive state given by Eqs. © with respect to small 
perturbation fields: 9(x,z,t), c(x,z,t), and 3>(x, z,t). 
The reduced set of three linear PDEs with constant co- 
efficients may be solved by the ansatz 

9(x, z, t) \ 

c(x,z,t) =uo(z) e lqx e at + c.c, (17) 



with the vector function 

/ \ 

u (z) = c(z) (18) 
V Hz)/(iq) J 

and c.c. denoting the complex conjugate, leading to a 
boundary eigenvalue problem with respect to z and the 
eigenvalue a. We solve the remaining linear ODEs by 
two different methods as summarized in the following 
paragraphs. 

The first one is the standard shooting method as de- 
scribed in detail for binary- fluid convection in Ref. [29| . 
The resulting coupled ODEs for the components of the 
vector function uq(z) are integrated for a set of initial 
conditions at one boundary. With the value of Uo at the 
opposite boundary a determinant f(a,R,q,Q,P,L,ip) 
follows. Keeping the initial conditions fixed, either R or 
a are varied such that / vanishes. The resulting values 
of a and R are functions of the remaining parameters. 

The second approach is based upon the so-called 
Galerkin method. The components of Uq(z) are expanded 
with respect to a suitable chosen set of functions fulfilling 
already the boundary conditions, i. e. either Eqs. (fT5|) 
or Eqs. (|16[) . Examples for this alternative numerical 
method may be found in Refs. (45l - |47| . The resulting 
generalized algebraic eigenvalue problem is solved numer- 
ically. 

At the onset of convection the small perturbations 9, c, 
and $ neither grow nor decay. This is the so called neu- 
tral stability condition, where the real part Re(cr) of the 
eigenvalue a with the largest real part (fastest growing 
mode) vanishes: 

Re(er)=0 with a = a{R,q,(,P,L 7 iP) . (19) 

This yields the Rayleigh number 

R (q) = R (q,C,P,L,^) (20) 

as a function of the chosen wave number q and describes 
the so-called neutral curve with a minimum at the crit- 
ical wave number q c and the critical Rayleigh number 
Rc = Ra(qc)- Convection sets in with a wave number 
q ~ q c by crossing R c from below. Depending on pa- 
rameters, the onset of convection may take place via a 
stationary bifurcation with a vanishing imaginary part of 
the eigenvalue, Im(cr) = 0, or via a Hopf bifurcation with 
a finite Hopf frequency Im(er) = ± uJo(q, C, P, L, ip) with 
its critical value oj c = wo (?<:)■ 

Throughout this work we choose for reasons of sim- 
plicity the Prandtl number P = 10. Since nanoparti- 
cles are much larger than, for instance, alcohol molecules 
in water, their mass diffusion is more than two orders 
of magnitude smaller. Accordingly, the Lewis number 
in nanofluids takes considerably smaller values of about 
L = 10~ 4 [48]. Growing linearly with the particle's size 
[48| the Soret effect can be changed in a wide range by 
varying the mass density with respect to the base fluid or 
the particle diameter. For small nanoparticles in water 
C ~ 0.01 and in Glycerin £ ~ 10 are appropriate values. 
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A. No-slip, impermeable boundary conditions 

It is a major result of this work that a particle- 
concentration-dependent thermal conductivity, described 
by finite values of (, leads in the range tp > to a shift 
of the onset of convection to larger values of R c . As an- 
other important result we find in the range tp ~ L and 
beyond some critical value £ c an exchange of instabilities 
from a stationary bifurcation to an oscillatory one and 
characterize this transition as a function of tp and £. 

The neutral curve Ro(q) belonging to the lowest bifur- 
cation from the heat conductive state at different values 
of £ is shown in Fig. [T] and Fig. [2] for two representa- 
tive positive values of the separation ratio tp, namely for 
tp = 1CT 5 « L (Fig. HJ and tp = 1(T 4 = L (Fig. ^. 
In both Figures the black solid line describes the neu- 
tral curve Ro(q) corresponding to molecular binary flu- 
ids, i. e., C — 0. It is included for illustrating the relative 
changes of the neutral curve as a function of £. We also 
mention for completeness that the Rayleigh number at 
the onset of convection in binary fluids decreases in the 
range of tp > with increasing values of tp starting at 
R c (tp — 0) ~ 1708. Simultaneously the critical wave 
number q c of the stationary bifurcation tends to zero in 
the range tp > L [H, H, H| . 

For increasing values of ( the neutral curve Ro(q) of 
the stationary bifurcation is shifted to larger values, as 
illustrated in Fig. [Ha). For C = 0.08 at about q ~ 3.55 
along Ro(q) a transition takes place from a stationary 
bifurcation in the range q < 3.55 to an oscillatory one 
in the range q > 3.55. We mention that here and in all 
following Figures solid lines mark stationary bifurcations 
whereas dashed, dotted, or dashed-dotted curves indicate 
oscillatory bifurcations. The neutral curve is determined 
by the condition given in Eq. (1191) for the eigenvalue 
with the largest real part. At about q ~ 3.55 the real 
parts of three different eigenvalues cross each other as a 
function of q, where the eigenvalue of the stationary bi- 
furcation has the largest eigenvalue in the range q < 3.55 
and a pair of complex conjugate eigenvalues leading to 
an oscillatory bifurcation has the largest real part in the 
range q > 3.55. Accordingly, at the transition between 
the two instabilities one has an unsteady change of the 
slope along the lowest neutral line. At a slightly larger 
value of £ this transition takes place at a smaller value of 
q and the path along the lowest parts of the two neutral 
curves belonging to the two instabilities becomes non- 
monotonic as illustrated in Fig. [T](b) for £ = 0.0875. At 
the crossing point of the two neutral curves in Fig. Hlb) 
the eigenvalue of the stationary branch vanishes as well 
as the real parts of the two complex conjugate eigenvalues 
at the oscillatory branch. This type of an exchange of in- 
stabilities is different from the well known codimension-2 
point in molecular binary fluids (2(| [2t| , where only two 
of the three eigenvalues are involved and one has a so- 
called double zero eigenvalue problem [3lj |. 

At larger, but still small values of the separation ratio 
tp the transition from a stationary to a Hopf bifurcation, 
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(a) Neutral curves for different f at tp = 10 -5 . 
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(b) Neutral curves near codimension-2 bifurcation. 

FIG. 1. Neutral curves Ro(q) are shown for different values 
of the conductivity parameter ( and tp = 10~ 5 . Solid lines 
mark stationary bifurcations and dotted/dashed lines oscilla- 
tory ones. In part (a) one has C — (black solid line), £ = 0.06 
(dark-gray solid line), £ = 0.08 (gray solid and dashed lines), 
and £ = 0.10 (light-gray dotted line). In part (b) two neu- 
tral curves, that correspond to the two eigenvalues with the 
largest real parts, are shown at £ = 0.0875. Again, the solid 
line indicates a stationary bifurcation and the dashed line an 
oscillatory one. 



which occurs by increasing £, takes place in a different 
and more diverse manner, as exemplified in Fig. [2] for 
tp = 10~ 4 . The neutral curve of the stationary bifurca- 
tion at ( = (black solid line) is strongly deformed by 
changing £ to £ = 0.05 (dark-gray curves). The mini- 
mum of the latter curve at q c is shifted to a smaller value 
compared to the black line. Moreover, for certain wave 
numbers, e.g. q = 2.5, one eigenvalue a vanishes at two 
different values of R by crossing two times the dark-gray 
solid line. Between these two values of R the growth rate 
Re(a) is positive, otherwise negative. Beyond the dark- 
gray solid curve at an even higher value for R at q = 2.5 
and £ = 0.05 the real parts of a pair of complex conju- 
gate eigenvalues change their sign and become positive, 
namely beyond the dark-gray dashed line. Increasing the 
conductivity parameter up to C = 0.15, the gray solid 
line for the stationary bifurcation results, which is even 
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(a) Neutral curves at ip = 10 4 . 

FIG. 2. Neutral curves Ro(q) are shown for tp = 10~ 4 and the 
heat conductivity parameter £ = (black solid line), £ = 0.05 
(dark-gray lines), £ = 0.15 (gray solid line), and £ = 0.20 
(light-gray dotted line). At £ = 0.05 the neutral curve splits 
into the dark-gray solid part for the stationary and the dark- 
gray dashed line for the oscillatory instability. At f = 0.15 the 
gray solid line marks the stationary branch and the oscillatory 
one is not shown. The light-gray dotted line corresponds to 
the oscillatory instability at £ = 0.20. 



further deformed compared to the case ( = 0.05. Increas- 
ing £ beyond C = 0.15 the stationary branch vanishes at 
about ( = 0.165 and beyond this value only a Hopf bi- 
furcation from the heat conductive state takes place as 
shown by the light-gray dotted line for £ = 0.2. Since 
the value of the Rayleigh number of the neutral curves 
at the Hopf branch changes only slightly as a function of 
C, the Hopf branch belonging to ( = 0.15 is not shown. 

This transition scenario, that increasing values of £ 
lead to a growth of the critical Rayleigh number at the 
onset of stationary convection and finally to a Hopf bifur- 
cation to convection at even higher Rayleigh numbers, is 
illustrated from a slightly different perspective by Fig. [3] 
In this Figure the critical values R c = Ro(q c ) 7 q c and uj c 
at the minimum of the lowest neutral curves are plotted 
as functions of the conductivity parameter £. The solid 
lines, describing the stationary bifurcation as a function 
of C, cease to exist at about £ = 0.165. During the pro- 
cess of disappearance of the stationary bifurcation the 
neutral curve of this bifurcation is deformed as illustrated 
in Fig. [2j Beyond C = 0.165 a Hopf bifurcation from the 
heat conductive state is preferred. Within the oscillatory 
region, R c varies only slightly with £ and takes values, 
that are more than 100% larger than the corresponding 
value for £ = 0. The critical wave number decreases 
initially, exhibits a major jump at the transition point 
and remains finally nearly constant, taking a value of 
q c ~ 3.2, which is close to the critical wave number for a 
simple Newtonian fluid, qf? F — 3.116. The critical Hopf 
frequency is monotonically increasing. 

A third perspective on the transition from a stationary 
bifurcation to a Hopf bifurcation is provided by Fig. 0] for 
£ = 0.1 and Fig.[5]for £ = 0.2. In both Figures we present 




(a) Critical Rayleigh number i? c (£). 




(b) Critical wave number q c {0- 




(c) Critical frequency w c (C)- 

FIG. 3. Part (a) shows the critical Rayleigh number R c , (b) 
the critical wave number q c , and (c) the critical frequency w c 
as a function of £ for ip = 10 -4 . Solid lines mark a stationary 
onset of convection and dashed lines an oscillatory one. 



the critical values R c resp. R c , q c and uj c as functions of 
the separation ratio i/>. As a guide to the eye we have also 
included the well known limiting case for £ = (gray solid 
lines). Again, solid curves describe the critical values for 
a stationary bifurcation, if this threshold is lower than 
the Hopf bifurcation, and dashed lines the critical values 
of the Hopf bifurcation. 

The left transition from a stationary bifurcation to a 
Hopf bifurcation in Fig. [4] and Fig. [5] is a codimension- 
2 bifurcation as illustrated in Fig. [TJb). The critical 
Rayleigh numbers at the minimum of the neutral curve 
of the stationary bifurcation at qf and at the minimum 
of the Hopf branch at qf? come rather close to each other 
with qf - ~ 0.025 in Fig. [5] 

With further increasing values of ip the transition from 
the Hopf branch back to a stationary branch is different 
for C = 0.1 (Fig. |U) and £ = 0.2 (Fig.©. The right tran- 
sition in Fig. 2] is similar to the left one and resembles 
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(a) Critical Rayleigh number R c {if>). 



FIG. 4. The critical Rayleigh number R c is given as a function 
of ip for £ = 0.0 (gray line) and £ = 0.1 (black lines). In the 
latter case the solid line marks the stationary and the dashed 
line the Hopf branch. At the transition from the stationary 
to the Hopf branch one has a codimension-2 point. 



qualitatively the scenario shown in Fig. HJb). Reducing 
£ below C = 0.1 the two transitions shown in Fig. H] ap- 
proach each other and one obtains a rather degenerate bi- 
furcation structure. The interesting nonlinear dynamics 
in this parameter range will be discussed in more detail 
elsewhere. The right transition in Fig. [5] shows a jump in 
the wave number and is in this respect rather different: 
By approaching this transition from larger values of ip 
the stationary branch ceases to exist (c.f. Fig. [2]). To 
the right of this second transition point convection sets 
in stationary. However, with decreasing values of ip the 
stationary part of the neutral curve becomes more and 
more deformed (just like in the case of increasing values 
of £ at a fixed ip, cf. Fig. [3]) and finally ceases to exist at 
the codimcnsion-2 point. Also this transition with quite 
disparate wave numbers at the two bifurcations bears an 
interesting nonlinear behavior to be discussed elsewhere. 

The ?/>-range where convection sets in via a Hopf bifur- 
cation changes as a function of the conductivity param- 
eter C as shown in Fig. [5] In the gray region within the 
displayed curve we find a Hopf bifurcation to convection, 
outside a stationary bifurcation takes place. Close to the 
left nose of this curve one has two transition points as in 
Fig. 31 whereas for larger values of £ the situation resem- 
bles the one shown in Fig. [5J with a strong jump in the 
wave number and the Rayleigh number along the upper 
part of the curve. 

In the limit of a vanishing Soret effect, i. e., ip — >■ 0, 
temperature gradients do not cause concentration gradi- 
ents anymore, the particle concentration becomes homo- 
geneous and the heat diffusivity \ is independent of the 
temperature. In molecular binary-fluid mixtures one has 
a codimension-2 point at i/'ctp — — L 2 , with a stationary 
bifurcation for ip > V'ctp and a Hopf bifurcation in the 
range — 1 < ip < ipcTP [ID, H(| . In colloidal suspen- 
sions with L ~ 1CP 4 or smaller this codimension-2 point 
V'Ctp is even closer to zero. 




(a) Critical Rayleigh number R c {iji). 




(b) Critical wave number q c (ip). 
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Ip 

(c) Critical Hopf frequency co c (ip). 

FIG. 5. In part (a) the critical Rayleigh number R c is shown 
as a function of ip for £ = (gray line) and £ = 0.2 (black 
lines), whereby in the latter case the solid lines mark a sta- 
tionary bifurcation and the dashed line a Hopf bifurcation. 
Part (b) shows the corresponding critical wave number q c {ip) 
and part (c) the frequency uu c (if}). 



In the range ip < ipcTP the Soret effect causes an en- 
hancement of the particle concentration near the lower 
plate and with 7 > also the heat conductivity in- 
creases at this boundary. For rising values of £ we find a 
monotonous reduction of the critical Rayleigh number at 
the Hopf bifurcation compared to the limit £ = 0. This 
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FIG. 6. Location of the two codimension-2-points (cf. Fig. [4] 
and Fig. [5]l in the (£-^)-plane. In the gray region inside the 
black curve the onset of convection takes place via a Hopf 
bifurcation and otherwise via a stationary bifurcation. 



is an opposite trend in comparison to the enhancement 
of the threshold found in the range tp > 0. However, 
this reduction of the Rayleigh number at the Hopf bi- 
furcation is rather small and no deformation of neutral 
curves takes place. With increasing £ the critical wave 
number q c decreases slightly and the frequency uj c in- 
creases almost linearly, similar to the behavior depicted 
in Fig. [3] (c). For example, at ip = — 10~ 4 and in the 
range < £ < 0.5 R c decreases with rising values of £ 
from R c = 1708 at £ = to R c = 1624 at £ = 0.5, q c 
decreases from q c = 3.116 to q c = 3.097 and uj c increases 
from u c = 0.19 to lo c = 0.28. 

The thermal conductivity contrast k is a function of 
the product RipC, = £. Since R c keeps large values at the 
Hopf bifurcation with decreasing ip, in contrast to the 
range tp > [26l [29J , the assumption of a linear variation 
of the particle concentration along the vertical z coordi- 
nate is violated at medium values of £, i. e., the linear 
dependence of the thermal conductivity on the particle 
concentration becomes violated for smaller values of £ 
than in the range ip > 0. 

The symmetry breaking caused by the nonlinear 
ground state of the temperature and concentration field 
generates characteristic changes of the flow field. Let us 
first concentrate on the case ip > 0. We found that for 
non vanishing values of £ the center zq of the convec- 
tion rolls is always shifted out of the center of the cell at 
2 = towards the lower plate, i.e., z$ < 0. For a fixed 
value of ip > the position zq decreases monotonically 
for rising values of £ as long as the onset of convection is 
stationary and, besides this shift, no further qualitative 
deformation of the convection rolls takes place. A typi- 
cal flow field for this situation is presented in Fig. Eta) 
for ip — 10~ 5 and £ = 0.085, where we show contour 
lines of the stream function $. This behavior changes, 
when £ takes even larger values for which convection sets 
in via a Hopf-bifurcation. In this paramter range we 
find inclined convection rolls, as exemplarily presented 



z 0.50 




(a) Stationary convection 




(b) Left travelling wave (ui c > 0) 



FIG. 7. Contour lines of the stream function <E> at the onset 
of convection for (a) tp ~ 10~ 5 and £ = 0.085 (stationary 
instability) and for (b) \& = 10~ 4 and £ = 0.3 (oscillatory 
instability). The center of the convection rolls is shifted to 
zo — —0.029 and zq = —0.005, respectively. Here, A = n/q c . 

in Fig. Efb) for tp = 10~ 5 and £ = 0.30. Left travelling 
waves (LTW, with tu c > 0) are inclined to the left, right 
travelling waves (RTW, with u) c < 0) are inclined to the 
right. Similarly inclined rolls appear also for negative 
values of the separation ratio and finite values of £, but 
in contrast to the results for tp > the center is slightly 
shifted towards the upper boundary for tp < 0. 

The dependence of the shift of the center of the con- 
vection rolls on the sign of the separation ratio tp can 
be understood as follows. For tp > 0, the suspended 
particles are driven to the upper plate of the convection 
cell and hence, the thermal conductivity rises from the 
lower to the upper plate. As the vertical heat current 
j z = — x<9 2 T con d is independent of z (compare Eq. [5} it 
follows straight forward, that the gradient of the tem- 
perature profile, which enhances the onset of convection, 
decreases from the lower to the upper plate. Therefore, 
the fluid motion tends to concentrate in the region of the 
lower boundary, which results into a shift of the center of 
the convection rolls to z$ < 0. Obviously, the situation 
is reversed for ip < and hence, the center of convection 
rolls is shifted upwards. 



B. Free-slip, permeable boundary conditions 

Let us now turn to free-slip, permeable boundary con- 
ditions. It is well known that for free-slip, permeable 
boundary conditions (|16p instead of no-slip, imperme- 
able boundary conditions (|15p the critical wave number 
q c does not tend to zero for increasing values of tp [26|, [2i| . 




300 I 1 1 ' ' 

1 1.5 2 2.5 3 

q 

(a) Neutral curves at ij) = 10 — 4 . 




12 3 4 

q 

(b) Neutral curves at ip = 10 — 3 . 

FIG. 8. Neutral curves Ro(q) for different values of the con- 
ductivity parameter £ are shown in part (a) at tp = 10~ 4 with 
C = (black solid line), C = 0.08 (dark-gray lines), C = 0.095 
(gray lines) and £ = 0.50 (light-gray dashed-dotted line) and 
in part (b) at ip = 10" 3 with < = (black solid line), < = 0.20 
(dark-gray lines), £ = 0.35 (gray lines) and £ = 0.50 (light- 
gray dotted line). Solid curves mark stationary instabilities 
and all ohter line types oscillatory ones. 

In spite of this difference the major trends are similar. 
For finite values of £ the threshold is enhanced compared 
to its value at £ = in the range tp > an d lowered 
in the range ip < 0. As in the case of no-slip, imperme- 
able boundary conditions we find a competition between 
a stationary and an oscillatory onset of convection in the 
range ip > 0, which, however, reveals some differences, as 
described in this section. 

Fig. [8] shows neutral curves for different values of £ at 
ip = 10~ 4 and ip = 10~ 3 . For smaller values of the sepa- 
ration ratio the major effect of finite values of £ is a shift 
towards higher Rayleigh numbers, very similar to the sit- 
uation depicted in Fig. [Ha). For ip = 10~ 4 the neutral 
curves are additionally deformed and split into a sta- 
tionary branch and an oscillatory one at higher Rayleigh 
numbers as illustrated in Fig. EJa). Even though the sit- 
uation is similar to the one shown in Fig. [2] there are 
some interesting differences. For £ = 0.08 the dark-gray 
solid curve merges with the dark-gray dotted curve with- 
out crossing, i. e. two stationary eigenvalues merge and 



9 



Re{a) 




350 400 450 500 



FIG. 9. The real parts of the two largest eigenvalues are shown 
as functions of the Rayleigh number R corresponding to the 
neutral curves in Fig. [8ja) for £ = 0.095 and at q — q c ~ 2.1. 
The black solid line marks the stationary branch and has two 
zero crossings, which dehne the boundary of the closed curve 
in Fig. [8ja). At about R ~ 502 this real eigenvalue merges 
with a second real eigenvalue (light-gray solid line) to a pair 
of complex conjugate eigenvalues (black dahed line). The 
appendant imaginary part is depicted in the inset. 



build a pair of complex conjugate ones similar as pre- 
sented in Fig. [9] for £ = 0.095. The island at finite values 
of q marked by the gray solid curve in Fig. EJa) and cor- 
responding to £ = 0.095 is a further feature which differs 
from the case of no-slip boundary conditions. The upper- 
boundary of this curve marks a linear restabilization of 
the heat conductive state, as indicated by the second sign 
change of the real part in Fig. [SJ After the restabiliza- 
tion an even higher positioned Hopf bifurcation occurs, 
as described by the gray dashed line in Fig. EJa) , which in 
the ranges of small and large values of q almost coincides 
with the oscillatory branch for £ = 0.08. As this restabi- 
lization takes place far beyond the first instability, where 
the nonlinear contributions to the equations of motion 
determine the dynamics of the system, an answer to the 
question whether this restabilization is of relevance for a 
real system can only be given by solving the nonlinear 
equations of motion in that range. 

For separation ratios larger than ip = 10~ 4 even a lin- 
ear restabilization of the oscillatory bifurcation may oc- 
cur. This is illustrated for ip — 10 -3 in Fig. |U[b). For 
£ = 0.2 we find a linear restabilization of the stationary 
branch, as depicted by the dark-gray solid line, before 
an oscillatory bifurcation from the heat conductive state 
(dark-gray dashed-dotted line) takes place. In contrast, 
for larger values of £, e. g., C = 0.35, the stationary 
branch has already disappeared and the lowest instabil- 
ity is given by a Hopf bifurcation. Here, this oscillatory 
branch shows a linear restabilization and a second Hopf 
bifurcation arises at even higher values of R as shown by 
the upper gray dashed curve in Fig. [Hfb). 

These interesting bifurcation scenarios including two 
exchanges of instabilities are presented from a different 
perspective in Fig. [TUJ Here, the critical values R c , q c , 
and uj c are shown as functions of £. In the first coexis- 
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(a) Critical Rayleigh number R C (Q. 
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(b) Critical wave number q c (C)- 
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(c) Critical frequency w c (C)- 

FIG. 10. The critical Rayleigh number R c in (a), the critical 
wave number q c in (b), and the critical frequency to c in (c) 
are shown as functions of £ for ip — 10~ 3 . Solid lines mark a 
stationary onset of convection and dashed lines an oscillatory 
one. 



tence range the situation is similar to the case of rigid 
boundary conditions, i. e., the end point of the solid 
line corresponds to the disappearance of an isola-like de- 
formed stationary neutral curve. Similarly, in the second 
coexistence range the first oscillatory instability ceases to 
exist at the end point of the dashed line, leaving behind 
a second oscillatory instability at even higher values of 
the Rayleigh number. We mention that unlike the end 
points of the solid and first dashed line the starting points 
of both dashed lines are not related to the creation or 
aniquilation of isolas. We decided to display these lines 
also to the left of the transition points to stress the fact, 
that the oscillatory branches are already present before 
another branch vanishes, cf. Fig. [3] The critical wave 
numbers and frequencies corresponding to the three dis- 
tinct instabilities differ considerably 

A further illustration of the bifurcation scenario is 



R c 




2.00 I — ■ 1 "— ' — ■ 1 "— ' — ■ — ' — • 1 "— ' — ■ 1 — J 

10~ 7 10~ 6 10~ 5 10~ 4 10 -3 , 

(b) Critical wave number q c (4>). 

FIG. 11. Part (a) shows the critical Rayleigh number R c as 
a function of the separation ratio ip for ( = (gray solid 
line) and £ = 0.2 (black lines), where in the latter case the 
solid lines mark the stationary branch and the dashed line 
the oscillatory one. Part (b) shows the corresponding critical 
wave number q c using the same line styles. 



given by Fig. 1111 where the critical Rayleigh number R c 
and the critical wave number q c are shown as functions 
of the separation ratio ip for a conductivity parameter of 
£ = 0.2. Again, we find the most pronounced changes in 
comparison to ( — in the range ip ~ L = 10~ 4 . When 
the two codimension-2 points marking the exchange of 
instabilities are plotted as functions of £, the resulting 
curve shows a similar behavior as the one for no-slip 
boundary conditions in Fig. [6] 

In accordance with the results for no-slip boundary 
conditions one obtains for negative values of the sepa- 
ration ratio ip and finite values of £ a reduction of the 
threshold compared to the case £ = 0. The critical 
wave number is independent from £ and takes the value 
q c = tt/\/2 just like in the case of molecular binary- fluid 
mixtures. 

For finite values of ( the changes in the flow field at 
the onset of convection with respect to the limiting case 
C = are qualitatively very similar to those discussed for 
no-slip boundary conditions. While a positive separation 
ratio ip > lowers the center of convection rolls, the 
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opposite happens in the case ip < 0. In addition, an 
oscillatory onset leads to inclined rolls. As for free-slip 
boundary conditions the velocity field takes finite values 
at both plates, the shift of the center is larger than in the 
case of no-slip boundary conditions. 



IV. DISCUSSION AND CONCLUSION 

The onset of thermal convection in a colloidal suspen- 
sion was investigated by means of a generalized contin- 
uum model for binary-fluid mixtures, which has been ex- 
tended beyond the Boussinesq approximation by taking 
into account a linear dependence of the thermal conduc- 
tivity on the local concentration of colloidal particles. 

We investigated colloidal suspensions where the ther- 
mal conductivity of the suspension increases with rising 
values of the particle concentration. An inhomogeneous 
particle concentration can be induced by temperature 
variations via the Soret effect, which corresponds to fi- 
nite values of the separation ratio ip. The concentration 
dependent thermal conductivity causes a nonlinear vari- 
ation of the temperature across the convection cell, in 
contrast to a linear variation for a constant thermal con- 
ductivity. The strength of the spatial variations of the 
heat conductivity is described by a dimensionless heat 
conductivity parameter £, as introduced in this work. It 
was found that for finite values of £ the vertical heat cur- 
rent through the convection cell is reduced compared to 
the case of a constant heat conductivity. 

For positive values of the separation ratio ip the sus- 
pended particles are driven to the colder upper plate of 
the convection cell and for a constant heat conductivity, 
i. e., £ = 0, a stationary bifurcation from the heat con- 
ductive state to convection takes place. Spatial variations 
of the heat conduction, corresponding to ( ^ 0, lead to a 
shift of the critical Rayleigh number R c to larger values 
than obtained in the case C — 0. Additionally, beyond 
some critical value £ c and in the parameter range ip ~ L 
the onset of convection takes place via a Hopf bifurcation. 
This trend, that a delay of the onset of convection leads 
to a Hopf bifurcation, is also met in the range ip < 0, and 
well known from molecular binary-fluid mixtures. 

Vertical variations of material properties, as discussed 
in this work for the thermal conductivity, are considered 
to be important for modeling convection in the Earth's 
mantle and, accordingly, a number of models including 
non-Boussinesq effects were explored 0, [l8| . For model- 
ing convection in systems with spatially varying material 
parameters also two superimposed layers of immiscible 
[49l - [5i| or even miscible fluids [52| are used. Like in our 
system with its spatially varying heat conductivity, an os- 
cillatory onset of convection was found in such two-layer 
systems for various parameters. 

The range of positive '(/'-values where the Hopf bifur- 
cation takes place increases with rising values of £. At 
both boundaries of this range we find a codimension-2 bi- 
furcation, where the thresholds of the stationary and the 



oscillatory bifurcation coincide. Additionally, one has a 
further codimension-2 bifurcation at ip < 0, which is well 
known from earlier investigations of molecular binary- 
fluid mixtures. However, for no-slip boundary conditions 
the two codimcnsion-2 points occurring for ip > are 
qualitatively different. Here, the two neutral curves be- 
longing to the two different bifurcations cross each other 
and therefore three eigenvalues are close to be critical. 
In contrast to this, only two eigenvalues are critical close 
to the codimension-2 bifurcation at ip < 0. At the left 
boundary of the ^-range, where a Hopf bifurcation is pre- 
ferred, one encounters usually small jumps in the critical 
wave number, whereas large jumps occur at the right 
boundary of this range. 

During our analysis we assumed 7 > in Eq. (fT]), cor- 
responding to the situation that the thermal conductivity 
of the particles is higher than the one of the base fluid. 
Hence, the local heat conductivity increases with a rising 
particle concentration. However, we obtain for 7 < ex- 
actly the same bifurcation scenarios as described in this 
work for 7 > 0. This may be understood as follows. In- 
dependent of the sign of 7, one obtains in the presence of 
thermophoresis layers of higher and lower thermal con- 
ductivity accompanied by a nonlinear z-dependence of 
the temperature. Changing the sign of 7 changes both, 
the variation of the heat conduction and the temperature 
profile, but the vertical heat current is identical for both 
signs of 7, as indicated by Eq. |5]). 

A further interesting question is how a spatially vary- 
ing heat conductivity affects the ratio between the con- 
vective and conductive heat transfer (Nusselt number) 
beyond the onset of convection, and how to understand 
the effects of boundary layers with enhanced or reduced 
heat conductivity. 
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Appendix: Determination of the heat conductive 
state 

In the heat conductive state, i. e. v = 0, the temper- 
ature distribution T coa d(z) and the particle distribution 
A r C ond(z) are determined by the two equations: 

= a,((l + 7(A r co„d - N ))d z T cond ) , (A.la) 

= d 2 z N cond + 3 2 z T cond . (A.lb) 
Jo 
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For a vanishing mass current at the boundaries, cf. 
Eq. (|4bf . a double integration of Eq. (|A.lb[) yields: 



termined by the two boundary conditions (|4a[) for T con d 
and one obtains: 



N cond = -—T cond + No + n . (A.2) 
Jo 



Now Eq. (lA.lap takes the form 

= d 2 z ({I + jno) T cond - ^rT 2 cond 



(A.3) 



T cond = r - J Ml ( i - ^ + Ml (\ + Z - ) , (A.5) 



2 d 



2 d 



with M± = T + T ± ST/2. The last free constant n in 
Eq. (|A.2[) is determined by the definition of Nq: 



and its integration gives 



1 



Ti ODd - T T cond = -(Ciz + Co) , (A.4) 

with To = (1 + "/no)To / ("/hx) ■ Solving this quadratic 
equation the two unknown co nstants Co and Ci are de- 



d/2 d/2 

N o = ^ J Ncond dz = N + n - J Tcond dz 

-d/2 -d/2 



(A.6) 
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